#!/bin/bash
set -e

function info() {
echo Usage: `basename $0` [-l in.bed] in.bam [in2.bam ...]
exit 65
}

while getopts  ":l:p:" opts
do
        case  $opts  in
        l) bed=$OPTARG;;
		p) out_prefix=$OPTARG;;
		*) info;;
        esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi


. /mnt/ilustre/app/medical/tools/.var



for item; do
input="$input -I $item"
done
echo $input

if test -n "$bed"; then

echo $bed;
echo;echo;echo gatk HP
java -Djava.awt.headless=true $tmp -Xmx$java_memory -jar $gatk \
-R $ref_genome \
-T HaplotypeCaller \
$input \
-o $out_prefix.vcf \
--dbsnp $data_path/ncbi/dbsnp/All_20150605.vcf.gz \
-L $bed
-stand_call_conf 30 \
-stand_emit_conf 10

else

echo no bed
echo;echo;echo gatk HP
java -Djava.awt.headless=true $tmp -Xmx$java_memory -jar $gatk \
-R $ref_genome \
-T HaplotypeCaller \
$input \
-o $out_prefix.vcf \
--dbsnp $data_path/ncbi/dbsnp/All_20150605.vcf.gz \
-stand_call_conf 30 \
-stand_emit_conf 10

fi


exit 0

	 
# Single-sample GVCF calling on DNAseq (for `-ERC GVCF` cohort analysis workflow)

java -jar GenomeAnalysisTK.jar \
     -R reference.fasta \
     -T HaplotypeCaller \
     -I sample1.bam \
     --emitRefConfidence GVCF \
     [--dbsnp dbSNP.vcf] \
     [-L targets.interval_list] \
     -o output.raw.snps.indels.g.vcf
	 
	 
# Variant-only calling on DNAseq
   java -jar GenomeAnalysisTK.jar \
     -R reference.fasta \
     -T HaplotypeCaller \
     -I sample1.bam [-I sample2.bam ...] \
     [--dbsnp dbSNP.vcf] \
     [-stand_call_conf 30] \
     [-stand_emit_conf 10] \
     [-L targets.interval_list] \
     -o output.raw.snps.indels.vcf
	 
	 
	 



Caveats
We have not yet fully tested the interaction between the GVCF-based calling or the multisample calling and the RNAseq-specific functionalities. Use those in combination at your own risk.
Many users have reported issues running HaplotypeCaller with the -nct argument, so we recommend using Queue to parallelize HaplotypeCaller instead of multithreading.